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Abstract 

A general reconstruction and calibration procedure for tracking and 
wire position determination of the OPERA drift tubes [T] is presented. 
The mathematics of the pattern recognition and the track fit are ex- 
plained. 



1 Introduction 



The two most important parameters for a track reconstruction with drift 
tubes are the time-distance-relation (rt-relation) to convert the measured 
times into spatial distances, and the resolution function as a weight for the 
track fit. In this article the determination of these parameters is described 
for the OPERA drift tubes. 

OPERA is a long-baseline neutrino oscillation experiment [2] to search 
for the Vn — > v T oscillation in the parameter region indicated by previous 
experiments [3J. The main goal is to find v T appearance by direct detection 
of the t from v T CC interactions. The OPERA detector consists of two 
massive lead-emulsion target sections followed by muon spectrometers. The 
task of the muon spectrometers, which consist of dipole magnets, RPC^] and 
~10000 aluminum drift tubes of 8 m length [1], is to clarify the signature of 
the muonic r decay and to remove the background originating from charmed 
particles produced in //-neutrino interactions. 

It is a difficult task to obtain the parameters needed for the track recon- 
struction of the muons traversing the OPERA drift tubes due to the very 
poor statistics of cosmic muons at the Gran Sasso laboratory where the 
OPERA detector is located. The determination of these parameters can, 
however, be done outside the underground lab using an equivalent detector 
with higher track statistics at same conditions (temperature, pressure, etc.). 
Since the mean cosmic energy at the sea level (~ 4 GeV) is the same mean 
energy of the muons coming out from the tau lepton decay after the CC 
v T interaction in the target, a calibration setup was built and operated at 
Hamburg. For this two 8 m long drift modules (pQ) with an total overlap in 
horizontal position were used, while the resulting wire sag was compensated 
by bending the tubes accordingly. This article describes a calibration proce- 
dure to determine the ri-relation and the resolution function. Furthermore 
a wire alignment procedure is described. 

2 The pattern recognition 

The pattern recognition fulfills two important functions. On the one hand 
it performs the preselection of hit candidates belonging to a track, on the 
other hand it delivers the start values for the track fit. 

We assume that a particle passing trough the detector (straight track) 
will produce TV" hits with valid drift times t within a defined time window. 
These times can be converted into distances rj(t) between the track and 
the fired wire by using the ri-relation. This radius ri(t) of a circle around 
the sense wire is an unsigned quantity and does not tell us on which side 
the particle passed the wire. The pattern recognition has to select the best 
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sample of hit candidates, define the signs to resolve the right-left ambiguities 
and define the start angle and the start distance to the origin of the track. 
For the calibration this procedure has to be fast and efficient. Normally a 
Kalman filteJU technique [2] is used. But if a sufficient^ number of events is 
available a more simple method described as follows fulfills the requirements 
above. It uses single straight tracks only. In case of more than one track 
per event (multi tracks) the event is ignored. 

For the track reconstruction with drift tubes the distance from the track 
to the wire is needed. Thus it is convenient to describe the particle track 
with the Hesse form 

do = x sin0 — y cos <f> . (1) 

Here do is the track distance of the closest approach to the origin and (j> 
the angle between track and x-axis as defined in Fig. [TJ In this analysis 
we assume parallel wires, therefore a two dimensional track description is 
sufficient. The distance of closest approach to the anode wires di is then 
calculated by 

di = do — Xi sin <fi + yi cos cj). (2) 

The index i describes the wire number i with the coordinates Xi and yi. 

In the simplest pattern recognition scheme the two tubes with the max- 
imum distance between each other in one event will be used (Fig. [2]). Now 
four possibilities exist to fit tangents (t± to t&) to the radii r± and ri- The 
tangent minimizing the \ 2 expression will be selected 

x =22— — ■ ( 3 ) 

2 A Kalman filter is a stochastic state estimator for dynamic systems. The most prob- 
able start value will be predicted and afterward corrected with the measured one. The 
difference of both values will be linearly weighted and is used to improve the current state. 

3 The sufficient number depends on the system. For the OPERA drift tubes 5000 tracks 
are enough for the calibration. 




Figure 1: Track description using the Hesse form and definition of the pa- 
rameters in the used coordinate system. 
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Figure 2: Example of the four tangents to the radii r\ and T2- The tangent, 
minimizing Eq. is the best description of the true track. 

a is the mean resolution, which is assumed to be the same for all tubes in 
the pattern recognition. For starting a is set to 1. The minimal x 2 °f t ne 




Figure 3: Sign convention for the track description considering the cases A 
and B according to the particle track on the right or left side from the wire. 

selected tangent has to be lower than a predefined maximal value (for this 
calibration values between 100 and 500 were used). For higher values the 
whole event will be rejected. The origin of such events can be noise or cross 
talk. If a valid tangent is found, the start parameter do start and (j) start are 
now available for the track fit. The signs of the di will be defined using the 
convention shown in Fig. El The normal vector ft is defined as 



:::J- w 




The particle can pass the wire at two sides, denoted by A and B. Eq. [I] can 
be written as 

d A B =™ B = Sil1 4> Xi A .B ~ COS ^ ViA B " ( 5 ) 
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Here %i A B and yi A B are the points of contact to the tangents. The distances 
di follow this convention, the vector x = Q') in Eq. [5] will be replaced by a 
vector pointing from the anode wire to the track 

3 Track fit 

3.1 General track fitting 

In this section a general track fit procedure will be described. It was origi- 
nally developed and used in the ARGUS experiment [5]. A former descrip- 
tion for the outer tracker system of the HERA-B experiment can be found 
in [9]. 

The distances d m ^ from the pattern recognition @ form a vector d m with 
the dimension N. The measurement uncertainties ai for each used wire are 
collected in the resolution function, binned in time or space respectively. The 
squared values of this function are on the main diagonal of the N x N co- 
variance matrix V. In the case of independent measurements only the main 
diagonal is filled. The track will be described by the N dimensional vector 
dt(q). In general q is M dimensional [M is the number of track parameters, 
in our case two). The parameter q results from the fit of dt(q) to d m . Using 
the least square method, the parameter q minimizes the expression 

d m -d t {q) W d m -d t {q) . (7) 



x 2 



The weight matrix W is the inverse covariance matrix of the measured 
coordinates: 

W = V -1 . 

Assuming independent measurements equation [7] can be written as: 

N 1 

^E*-^!))- ( 8 ) 

i=l i 

Minimizing \ 2 as a function of q is equivalent to finding a solution to the 
system of equations 

% = °- < 9 > 

In general these equations are non-linear because the track model d{q) 
is non-linear. However, in most of the cases it can be linearized and the 
equations can be solved iteratively. A linear track model is given by: 

d t {q) = Aq + b with A = (A itl ) i=i,..., N . (10) 

fj.=\,...,M 



The index m denotes the distances calculated from the measured drift times. 
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Here A denotes the Jakobi matrix. With this, Eq. Ogives 



[d m - Aq-Pj w(d m -Aq-b 



X 



(11) 



N 



M 



—{dm,i — X! Aikqk — h) 



i=i 1 



By differentiation one gets 

dx 2 
dqk 

dx 2 
dq 

This yields the parameters 



N 



k=l 



1 



2 ^ ] A-i^ 2 {dm,i ^ ] Aikqk hi 
fc=l 

A?- 6) 



i=l 

-2 [A T W (d 



-2A T W (cT m - Ag - 5) D = J 0. 



(12) 
(13) 



(A T WA) 1 A T W (d m - Pj . 



(14) 



With the found parameters one can solve the following equation system: 



A T WAq = A T W [dm - Pj . 
From the general formula for error propagation ([B]) 

V(y) = ZV(x)Z T for y = Zx + c 

and the replacements 

x = d m 
V(x) = W" 1 

y = ? 

Z = (A T WA) _1 A T W 

one receives the covariance matrix 



(15) 



(16) 



V(<f) = (A 1 WA)- 1 . 



(17) 



As already mentioned this method can be used for a non-linear track model 
if it can be linearized. A linearisation requires a starting value % near the 
true solution. In this case one can describe the vector %, which essentially 
contains the non-linearities, with the first order of the Taylor expansion 



M 



(qh - = dt(oo) + Mq- <fo) 



(18) 



q=q a 
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The Jakobi matrix of dt(q) for q = % is defined as 



ddt 



(^) 



<2=<?o 



ddtA 

^Qft J i=l N 

H=1,...,M 



(19) 



Correspondingly b = dt(qo) — A% (see Eq. flOl) can be used in Eq. [15] 
This results 



(20) 



(21) 



A T WA(?- <f ) = A T W [d m - d t (Qo) ■ 
The change of parameters 

Aq = q- q 

results into a better approximation for the next iteration 

q\ = q + Aq. (22) 
Eq. [20] and [22] define an iterative procedure. In the nth iteration, equation 

A T WA(Ag) = A T W \d m - d t (q n ^)} . I (23) 



is solved for Aq. The vector d m {q n -i) contains the fitted track distances for 
the track parameters of the next iteration and A is the Jakobi matrix of the 
track parameter for q = q n -\. Aq then gives the new approximation of the 
parameters 



q n = q n _ 1 + Aq. 
The end condition is reached when Aq ~ and 



An 



Xn-1- 



(24) 



(25) 



3.2 Track fit with two parameters 

In this section a minimal track fit with two parameters will be described in 
detail for the drift tube modules. The initial parameter qstart for the track 
fit found by the pattern recognition is 



'd 



qo = Qstart 



Ostart 



(26) 



tart 



The Jakobi matrix (Eq. [T9]) can be written as: 



dd 



q=q n -i 



ddA 

dq^ ) i=i,..., N 

1 |U=1,2 



(27) 



7 



{Aifj) contains the partial derivatives of di (see Eq. [2]) with respect to do 
(fi = 1) and <f> (fj, = 2) 

ddi ddi 

—j- = 1 and — = -Xi cos0 - yi sm<p. (28) 
ado ocp 

By substitution (Eq. 1251) 

£? = A T WA and y = A T W(c? m - t£(^_i)), (29) 
one receives 

GAq = y and = g~ x y. (30) 

5 is an 2 x 2 matrix while y and g are two dimensional vectors. For inde- 
pendent measurements Q and y can be calculated by 

^ 1 ddi ddi 1 , Q1 s 

^ = ^^ = 1 ' 2 ' (31) 



^ = E^(<W-<M- (32) 
i=i a ^ CT « 

The calculation of the partial derivatives is done in Eq. [271 The subscripts 
ji and v describe the elements of the parameter vector q. Calculating Aq 
with Eq. [30] gives the new iterated parameters 

q n = q n -i + Aq. (33) 

This procedure (Eq. IHU|) will be repeated, until the end condition 

Ax 2 = xl-i ~xl< AxLn (34) 

is fulfilled. A quality check for the track fit is the distribution of the x 2 
probability, shown in Fig. [H On the left side events have been clustered 
with an under-estimated resolution whereas the resolution has been over- 
estimated for the events on the right side. Except for the peak on the left 
side due to noise and cross talk, this distribution is flat, since a uniform error 
distribution is expected. Furthermore, the mean value of the x 2 distribution 
per degree of freedom, also shown in Fig. [H should be around one if the 
track fit works correctly. In reality the mean value varies between 1 and 1.3. 

For the fit described here a Ax^ in = 10~ 7 was used. The procedure 
converges usually after three iterations. After that the signs of d m> i, de- 
termined by the pattern recognition, will be corrected (turned) if they are 
opposite to the signs of dt i . This procedure is repeated until no sign changes 
appear anymore (solution of left/right ambiguity). Such hits with random 
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Figure 4: Probability distribution of \ 2 (left) and the distribution of x 2 P er 
degree of freedom (right). 

time signals, coming from noise or cross talk, are removed if the \ 2 distri- 
bution is too large. The whole procedure is repeated, starting with the qo 
from the pattern recognition, until \ 2 < Xmax- After this it will be tested if 
still enough tubes N are left over for the reconstruction (> N tu bes min = 4). 
In this case, the vector q n of the last iteration contains the track parameter 
with the best track description. If not, the whole event is rejected. 

4 Calibration 

The calibration is the iterative determination of parameters like the rt re- 
lation and the resolution function a (to)- In addition the alignment of the 
wire positions is possible with this procedure and will be described. 

The calibration procedure is based on the concept of residuals e(io), 
defined as difference of the fitted drift distance (dt) and the measured one 

(d m ) 

e(t D ) =d t - d m . (35) 

The sign of d m is defined as sign of dt- The rt-relation and the resolution 
function are represented by subdividing the parameters in time bins for an 
iterative procedure. 

4.1 Time distance relation 

The correlation of the drift time to and the drift distance is given by the 
rt-relation using the drift velocity vd 

d t (t D )= / v D (t)dt= / —dt. (36) 
JO Jo at 
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Figure 5: The signs of the drift distances, the residuals and the corrections 
of the rt-relation. 

The iterative correction of the rf-relation is done by the difference Ad m (tu), 
calculated from the residuals including the signs of the drift distances 

Ad m (t d) = signidt) e(i D ) = signid t )id t - d m ). (37) 

Admits) specifies, whether the measured drift distance is too low or too 
high. The signs used in the above equations are presented in Fig. The 
drift distance dt(io) is defined negative, if the particle track passed the 
wire on the left side (according to the defined coordinates) and positive on 
the right side. The sign of the residuals, arising from the drift distances, 
follows the same convention if the distance of the fitted track is larger than 
the measured drift distance. In the other case, the sign of the residuals is 
opposite to the drift distances. If Ad m (io) > 0, the rt-relation is increased 
in the given time bin for the next iteration, for Admit d) < it is decreased. 
Normally the residuals are Gaussian distributed. In the vicinity of the wire 
the influence of the primary ionisation statistics increases due to the cluster 
distribution, and the drift time distribution will be asymmetric to larger 
values and thus to larger drift distances. Using Eq. [37J the distribution of 
Admit d) will be asymmetric to the negative side. This effect will be stronger 
the closer the track passes the wire. Since the difference in Eq. [35] can not 
be smaller than the negative value of d m (same sign as dt), this effect leads 
to a cut on the negative side in the Admit d) distribution (see chapter |3~2|) . 
The cut on Admit d) distribution is exactly — \d m \. Thus the distribution is 
positive at the wire position. In all cases the maximal entries are at zero 
position. The mean value of Admit d) corrects the rt-relation as follows 

d m itD)new = d m itD)old + [Admit d)), (38) 

with [Adm[t d)) as mean value of the Ad m (io) distribution. The iteration 
stops if the squared sum over all Admit d) is n °t changing with further 
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Figure 6: Correlation between drift times and distances, called ri-relation 
(left). On the right side the distribution of the fitted distances are shown. 

iterations. 

For the start, the distance d m (tD)oid ls obtained by integration of the 
drift time spectrum. In a first approximation a drift tube with a homo- 
geneous distribution of track distances will be assumed. With dN tracks, 
passing the wire within the interval [r; r + dr], 

dN N track 

—— = const = , (69) 

dr r tube 

Ntrack is the number of tracks which hit the tube and r tu b e is the tube radius. 
The drift velocity vd is given by 

VD = rtube . — (40) 

Ntrack dt 

By integration of the drift time spectrum dN/dt and normalizing to the 
maximal tube radius one gets a first estimation of the rt-relation 

, x r tube f tD dN 

d m ,start{tD) = -T- T / ~~T~ dt. (41) 

Ntrack JO at 

Fig. [6] shows the correlations between the distances dt{to) and the drift 
times tp. The mean values of the dt(to) distributions of both branches 
form the ri-relation. For the calibration two versions of rt-relations are 
used. The description of the time binned version is a good approach for 
distances not too close to the wire. Since the residua distribution near the 
wire is asymmetric as mentioned above, the rt-relation (time binned) is 
bend away from the actual wire position. This description is unrealistic, 
since the position zero is existing. To avoid this behavior, a combination of 
time binned and space binned ri-relation is used, where the space binned 
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relation is used for positions near the wire. Thus the wire position is reached, 
but the drift times are always greater than zero. This behavior is realistic, 
since a drift time zero can not appear due to the cluster distribution even 
in the case if the wire is hit by a track. For describing, interpolating and 
smoothing the rt-relation cubic splines are used. If the method is working 
correctly, the distribution of the fitted distances is theoretically flat. In 
reality this distribution is a bit bended, but should be without spikes and 
gaps. 

4.2 Spatial resolution 

For the determination of the resolution function ajitn) the RMS value of 
the residual distribution eiitn) will be used 



N specifies the entries of each residual distribution for a given time bin j and 
fi is the residual scale factor. Since the residuals are calculated as the differ- 
ence of measured and fitted track distances, they contain the uncertainties 
of the fitted track. To avoid underestimating the resolutiorS a scale factor 
for correction is needed. The calculation of this factor will be described in 
the following [7J. 

The resolution function is used to construct the covariance matrix V((f m ) 
which appears in the x 2 expression (Eq. [7]). The expectation of the squared 
deviations of measured and fitted distances to the wires are the diagonal 
matrix elements 



Since per definition the true track (parameter qt rue ) has no errors, the co- 
variance matrix can be written as 



That means the covariance matrix elements are equal to the standard de- 
viation of the true residual distribution. Since the true track (and the true 
residual distribution) is unknown, it will be replaced by the parameter q of 
the fitted track of the last iteration. The relation of the measured residuals 
to the covariance matrix is 




(42) 




(43) 



V(d m ) = V(d m - d t {qtrue)) = V(e* true ). 



(44) 



V(e) = V(d m - d t (q)) + V(d) 



(45) 



By linearisation 




(46) 
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Only the fitted track is known, not the real one. 
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the true residuals can be expressed by measurable values 



d t {q)+A(q true -q) = e - A(q true - q) . (47) 



(■true — d m dt((lt rue ) — d n 

By means of a general error propagation (|6j) one gets 

V(el rue ) = V(e) + V(Aq) = V(e) + AV(<f)A T . 



(48) 



This means, the reference trajectory dt(q) is only known within the track 
parameter errors. If only the standard deviation of the measured residual 
distribution is used for the correction of the resolution function, it will be 
underestimated (further information in [H]). A convenient correction is done 
by the residual scale factor. With the requirement 



Vu&rue) D = f f! ■ V^e) = /? 



Vn(et r 



(AV(<f)A 5 



one gets the factor 



fi 



at 



(49) 



(50) 



for the ith measured residual ej of the track. A is the Jakobi matrix (Eq. 
|27|) and V = W _1 the covariance matrix of the track parameter vector 
q (Eq. |50|) . The field cr, contains the calculated resolution (Eq. [4*2]) for 
the given drift time bin of the last iteration. Fig. [8] shows the residual 
distribution, corrected with the scale factor. If the complete covariance 
matrix is unknown, a global correction factor can be used. Assuming that 
all N fired tubes have the same Gaussian distributed error a and the same 
mean squared measured residual e 2 , one gets 

■2 



N- 



(7 



ndof 



N 



(51) 



with ndof degrees of freedom (2 for a linear track model). 

Each residual €i(tn) will be corrected by the factor and filled into his- 
tograms. The RMS value of this distribution is the mean resolution a{t£>). 
Repetition for each time bin i yields the resolution function ai(tn) depen- 
dent on the drift time. The weight requires a resolution for each individual 
wire distance d and d±. For the interpolation and smoothing cubic splines 
are used. Thus resolution values for a finer binning are available. Fig. [7] 
shows the resolution function and the splines. The time dependent resolu- 
tion function can be converted into a distance dependent one by using the 
rt-relation. 

A quality factor of drift tubes is the mean resolution 

(") = , = , (52) 
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Figure 7: Resolution function (squares) dependent on the drift time, fitted 
with cubic splines (diamonds). 

Reco plots for all events, run data_031_08-55-13.operaasc 
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Figure 8: Distribution of the residuals corrected by the scale factor. The 
RMS value yields the mean spatial resolution. 
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with <7j as resolution of the ith of N drift time bins (e.g. 20 ns per bin). 
This definition takes into account, that distances with a better resolution 
have a higher weight in the track fit. Using the mean resolution (a) in the 
next iteration of the track fit, the mean errors of the track parameters are 
equal to the weighting of a variable resolution function a(tn) [Z]- 

The resulting resolution function can be described by a theoretical model 



It consists of contributions from the primary ionisation statistics, from the 
diffusion and from a constant part k (uncertainties in wire position, time 
jitter in the electronics, time resolution, etc.). The diffusion, dependent on 
the field and the drift distance, is specified by d. The parameter j is the 
number of the triggered cluster, np the number of primary ionisations per 
unit and r the drift distance. 

4.3 Wire calibration 

The presented procedure can also be used for the calibration of the wire po- 
sitions by replacing the track parameter vector q with a vector (3 describing 
the wire positions 




(53) 




(54) 



The track description is then given by 



di = d - fix sin(0) + [3 y cos (<f>). 



(55) 



The Jakobi matrix can be written as 




(56) 



with 



— — = -5 lj sin <j) and — — = <% cos 



(57) 




i and j denote the wire number. Equation [231 can be written as 



A T WA(A/3) = A T W \d m - d t 0n-i) 



(58) 



A/3 contains the wire shifting in the x and y direction. The calculation of 
A/3 is done by 

GAp = y and A/3 = ^" 1 ^, (59) 
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with the determination of Q and y as in Equations [29J [31] and [32] by replacing 
q with (3. The stop condition is the same as described by Eq. [M] In this 
case A/3 should not change (after 2-3 iterations). For a start parameter 
the theoretical wire positions are used. 

In most cases, e.g. for a narrow angle distribution, it is more sufficient to 
align only one wire coordinate (perpendicular to the track axis). A simple 
way is an iterative procedure using the mean value of the residual distribu- 
tion. For each iteration the wire position will be corrected by the deviation 
of this mean value from zero, attenuated by a scale factor of 0.3 to 0.5. 
It is of advantage that not only tracks perpendicular to the wire deviation 
have to be used. Tracks at angles around those also give good results and 
act as additional attenuation. This method is very stable, but without a 
scale factor the technique is not working and the wire correction starts to 
oscillate. Furthermore this method only works for homogeneously irradi- 
ated drift tubes, for tubes at the edges the position correction is biased and 
should not be used. Fig. [9] shows some plots of the wire positions before 
(upper plots) and after the correction (lower plots). In the left plots the 
deviation of the wire position is shown, in the right plots the distribution 
of the deviation can be seen. The RMS value gives the overall mean wire 
misalignment. The plotted results show the improvement of the accuracy of 
the wire position from 183 /im down to 38 /mi after five iterations. 

4.4 Calibration schematic 

A flow diagram of the iterative calibration procedure to determine the rt— 
relation and the resolution function is presented in Fig. [101 For the first 
iteration a rt-relation from the integration of the TDC spectrum is used. 
Furthermore, the run time correction (to), not described in this article, can 
be extracted from the TDC spectra for each channel and used for the drift 
time calculations. For the first iteration a constant resolution function (a) = 
1000 /im is used in the pattern recognition and the track fit. The procedures 
of the pattern recognition and the track fit are described in detail in chapter[2] 
and[3] After the iteration the resulting rt-relation is used as start parameter 
for the next iteration. This procedure will be repeated until the rt-relation 
will not change. The indicator is the sum of the squared differences of old 
and new rt-relation. Now the procedure described above will be repeated 
using the resolution function. The convergation of the rf-relation and the 
resolution function appears after several iterations (typically 10). At the 
end a wire calibration can be done if necessary. In this case, the whole 
calibration has to be repeated, starting from the drift distance calculation 
(see Fig. [T0j) . After several iterations (typically 5), the wires are aligned. 
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Figure 9: Wire position before (upper plots) and after the wire alignment 
(lower plots). In the left plots the deviation of the wire position is shown, 
in the right plots the distribution of the deviation. 
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Figure 10: Schematic of the calibration. For the first iterations (dashed 
lines) a mean resolution is used. The wire calibration is done once for the 
determination of the wire positions. 
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